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Abstract 

Supernova explosions provide the most sensitive probes of neutrino propagation, such as 
the possibility that neutrino velocities might be affected by the foamy structure of space-time 
thought to be generated by quantum-gravitational (QG) effects. Recent two-dimensional 
simulations of the neutrino emissions from core-collapse supernovae suggest that they might 
exhibit variations in time on the scale of a few milliseconds. We analyze simulations of such 
neutrino emissions using a wavelet technique, and consider the limits that might be set on 
a linear or quadratic violation of Lorentz invariance in the group velocities of neutrinos of 
different energies, v/c = [1 ± (E / M v \x{)] or [1 ± (E / M v yw2) 2 ], if variations on such short 
time scales were to be observed, where the mass scales M„LVi might appear in models of 
quantum gravity. We find prospective sensitivities to M^lvi ~ 2 x 10 13 GeV and M uL y 2 ~ 
10 6 GeV at the 95 % confidence level, up to two orders of magnitude beyond estimates made 
using previous one-dimensional simulations of core-collapse supernovae. We also analyze the 
prospective sensitivities to scenarios in which the propagation times of neutrinos of fixed 
energies are subject to stochastic fluctuations. 
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1 Introduction 



Supernovae provide some of the most sensitive probes of neutrino physics pQ , as exemplified 
by studies of the neutrinos detected after emission from SN 1987a [2]. Excellent prospects 
for improving these studies would be offered by a future galactic core-collapse supernova. 
Crucial inputs into estimates of the prospective sensitivities of such studies are provided by 
supernova simulations. In the past, these simulations have mainly been one-dimensional, 
i.e., implicitly assuming that the collapse is spherically symmetric. Clearly, more realistic 
simulations are desirable, and more recently two-dimensional hydrodynamic models, i.e., im- 
plicitly assuming only cylindrical symmetry, have become available, in which multi-group, 
three-flavor neutrino transport has been treated by different (relatively sophisticated) ap- 
proximations [3H5]. These have revealed several interesting features, the most relevant for 
our analysis being the appearance of fast time variations in the neutrino emission [4j|5], as 
seen in the top panel of Fig. [TJ It was shown that these could be observable in the IceCube 
experiment [6], and could exhibit a quasi-periodicity of O(10) ms, similar to the natural time 
scale of reverberations associated with hydrodynamic instabilities in the supernova core. It 
is desirable to confirm the predicted appearance of such features in neutrino emissions from 
core-collapse supernovae through more detailed simulations, in particular in three dimen- 
sions 0. Nevertheless, we find these features sufficiently interesting and well motivated to 
consider the sensitivity to effects in neutrino propagation that would become available if 
such rapid time variations were in fact to be observed. 

The most obvious such effect is that of neutrino mass. Non-zero masses cause neutrinos to 
travel at less than the speed of light, by an amount that decreases with increasing neutrino 
energy. This causes any time structure that appears simultaneously in emissions over a 
range of energies to spread out before arrival at the Earth, an effect that was exploited to 
set an upper limit on neutrino masses using data from SN 1987a [9]. Though an interesting 
demonstration of principle, that limit was not competitive with laboratory and cosmological 
limits, and even the increase in sensitivity suggested by the more recent two-dimensional 
simulations seems unlikely to be competitive. 

Another possibility is to use the neutrino emissions from core-collapse supernovae to 
constrain effects on neutrino propagation such as might be induced by 'foamy' quantum- 
gravitational fluctuations in the fabric of space-time [10] • Such space-time foam effects could 
include a refractive index (i.e., a change in the neutrino or photon velocity that depends on 
energy) [TTHT5] . dispersion in propagation at fixed energy [20], and a loss of coherence pjjj. 
Models suggest that any such effects should increase with increasing neutrino energy E, being 
proportional, e.g., to E/M uL y 1 or (E / M uL y 2 ) 2 where the mass scales Mj, L vi might originate 
from quantum gravity (QG). As such, they would be easily distinguishable, in principle, 
from the effect on neutrino propagation of a neutrino mass. Estimates have been given [15] 
of the possible sensitivity to such effects that could be provided by a core-collapse supernova 
explosion in our Galaxy, e.g., sensitivity to M uL yi ~ 2 x 10 11 GeV or M uL y 2 ~ 2 x 10 5 GeV 
for a refractive index > 1, corresponding to increasingly sub luminal propagation of energetic 
neutrinos. These prospective sensitivities are considerably greater than those offered by 

1 For first steps in this direction, see pj[8]. 
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terrestrial long-baseline neutrino experiments [16J: the latters' beams have much finer time 
structures, but they are handicapped by their much shorter propagation distances. 

The estimates above were based on the earlier simulations of the core collapse of a 
supernova that yielded emissions over a period of seconds without any finer time structure. A 
priori, the observation of time structures on the scale of O(10) ms in the neutrino emissions 
from core-collapse supernovae, as suggested by more recent two-dimensional simulations |^]|5] , 
would provide the possibility to constrain foamy effects on neutrino propagation a couple 
of orders of magnitude more sensitively than was previously estimated. In this paper, we 
use such a simulation to estimate the prospective sensitivity to a neutrino refractive index 
and dispersion in propagation at fixed energy, applying a wavelet analysis to the simulated 
neutrino signals published in [5]. 

We find that the prospective sensitivities to novel effects in neutrino propagation would 
indeed be enhanced by two orders of magnitude if the neutrino signals from core-collapse su- 
pernovae do exhibit the fine-scale time structures suggested by two-dimensional simulations. 
In the first instance, these sensitivities may be expressed for energy- dependent time shifts r n 
or dispersions at fixed energy, both expressible in units of s/MeV' for effects oc E l . Knowing 
the distance from any given supernova, these sensitivities may be translated into sensitivities 
to model parameters such as the M^lvz introduced above. Hypothesizing a typical distance 
of 10 kpc we find, for example, sensitivities to M„lvi ~ 2 x 10 13 GeV and Mj,lv2 ~ 10 6 GeV 
at the 95 % confidence level. We emphasize these sensitivities could immediately become 
lower limits if time structures are observed. The absence of such time structures could pro- 
vide upper limits if the predictions by two-dimensional simulations of short time structures 
could be validated, assuming that no astrophysical effect during propagation could wash out 
the effect. 

This subject has become much more topical during the finalization of this paper, with the 
publication of an analysis of OPERA data [T7] using a technique similar to that developed 
in [15]. This reports possible superluminal neutrino propagation with velocity exceeding that 
of light, c, by an amount Sv: 5v/c ~ 2.5 x 10 -5 , corresponding to M^lvi ~ 1-1 x 10 6 GeV or 
^i^LV2 ~ 5.6 x 10 3 GeV [18] . smaller than the lower limits established in [15] using SN1987a 
data. Thus, a simple power-law fit jointly to the SN1987a data and OPERA Sv would prob- 
ably require an energy dependence steeper that E 2 , in apparent conflict with the energy 
spectrum of the neutrino events measured by OPERA. The superluminal neutrino interpre- 
tation of the OPERA data is subject to many other experimental and phenomenological 
constraints, and one should not assume that it will survive further scrutiny. Nevertheless, 
this episode heightens awareness of the importance of probing fundamental principles such 
as the universality of the velocity of light as sensitively as possible, and our analysis based 
on two-dimensional simulations of supernova explosions shows that they could provide un- 
parallelled sensitivity to novel effects in neutrino propagation. 
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2 Ingredients in the Analysis 



2.1 Quantum-Gravity Models for Non-Standard Neutrino Propa- 



As mentioned in the Introduction, various possibilities for non-standard neutrino propagation 
are suggested by phenomenological models based on approaches to quantum gravity [T5|[T9]. 
This is because such models entail microscopic fluctuations of space-time, due to curvature 
fluctuations and/or, in certain theories, space-time defects. Specifically, in brane theories 
based on string theory [TJ], the latter may be modelled as point-like structures that cross the 
brane Universe from the bulk, giving space-time a "foamy" nature at microscopic scales. In 
our current state of knowledge of string theory, the string scale is essentially a free parameter 
to be constrained by low-energy phenomenology, and in particular by constraints on non- 
standard neutrino propagation. 

2.1.1 Modified Dispersion Relation 

One of the most explored avenues for experimental probes of (some models of) space-time 
foam is to search for Lorentz violation induced in the propagation of matter particles by their 
interactions with this foamy space-time medium. In certain string-inspired models of foam 
the presence of a medium affects the dispersion relations of certain species of matter parti- 
cles |14[I20]. In the simplest formulation of the effects of the foam, in a first approximation 
only electrically-neutral particles interact with the medium, so that photons and neutrinos 
are the most sensitive probes of such models. In these models, the modification of a parti- 
cle's dispersion relation is a consequence of the microscopic Lorentz violation induced by the 
recoil of a space-time defect during its non-trivial interaction with the open-string state that 
represents the particle excitation in a brane Universe [HIEO]. In the case of photon propa- 
gation, this effect is manifested as a vacuum refractive index. For purely string-theoretical 
reasons, the induced refractive index is subluminal, implying that, if a beam of photons with 
different energies is emitted simultaneously from a source, the arrival times of more energetic 
photons will be delayed compared to their lower-energy counterparts @. 

In the case of neutrinos, regarded as (almost) massless particles, a similar effect might 
be expected, namely a delayed arrival of the more energetic neutrinos from cosmic sources, 
assuming that the neutrinos of different energies are emitted (almost) simultaneously. One 
may therefore consider foam-induced Lorentz violation that is expressed via a neutrino group 
velocity, v g , that may depend either linearly or quadratically on the energy of the neutrino: 



As discussed in [15] . lower limits on M„lvz m &y be obtained by requiring that narrow peaks 
in neutrino emission over a range of energies not be broadened significantly, or even washed 
out, and this is the strategy advocated here using emissions from a supernova. 

2 As a bonus, this avoids any constraints due to the emission of Cerenkov radiation by energetic pho- 
tons [21] . 
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2.1.2 Dispersion of a Wave Packet 

A second effect that may also be induced by quantum-gravity foam models is a spread in 
the width of the wave packet, which may depend on the neutrino energy. This could arise 
from an energy-dependent neutrino velocity of the type discussed above, or from a stochastic 
spread in neutrino velocities at fixed energy [20] . 

Consider, for example, a neutrino wave packet that is Gaussian in the "approximately" 
light-cone variable x — v g t, where the group velocity v g = du/dn is near the speed of light 
for the relativistic, almost massless neutrinos that we consider here, allowing for a generic 
dispersion relation u = with k = \k\ being the spatial momentum amplitude. For 

neutrinos that are almost massless, the analysis is similar to that given in [20]. The square 
of the neutrino amplitude is given in general by: 

l/0M)| a = ; ^ = exp I (X ~/f a2t2 i > , ( 2 ) 



(AxoV 



where \Axq\ is the spread at t — 0, a — ^ {d 2 ui /cPk), and we assume that the neutrino wave 
packet has a Gaussian distribution in the "approximately" light-cone variable x — v g t. 

We see immediately in ([2]) that the quadratic term a in the dispersion relation does not 
affect the motion of the peak, but only the spread of the Gaussian wave packet: 



a 2 t 2 

+ (3) 

which therefore increases with time. The quadratic term a also affects the peak amplitude 
of the wave packet: the latter decreases as the spread ([3]) increases, in such a way that the 
integral of \f(x,t)\ 2 is constant. 

If the neutrinos were exactly massless, then, as in the case of photons, the quantity a 
would receive non-zero contributions only from anomalous terms in the dispersion relation 
due to quantum gravity of the form p]) [20]. In the linear case, I = 1, such corrections would 
be independent of the energy of the neutrino. In the presence of small neutrino masses, 
m<K, there are always contributions to a from terms of the form: 

2 

m 

a = — , m <C k . (4) 

K 

Such terms contribute to the spread of the wave packet, but decrease with the neutrino 
momentum, in contrast to quantum-gravity effects that are expected to be constant or to 
increase as functions of the neutrino energy (momentum), depending whether the energy 
dependence of the refractive index is linear or quadratic. 

Hence we may parametrize the spread a generically as in fl3]), with the parameter a 
having a power-law dependence on the neutrino momentum k: 

i i— i 

a "~ /(Z + 1) M^' (5) 
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with I = 1 or 2, where the stochastic case is denoted by the tilde in the suffix of the QG scale. 
To leading order in the small neutrino mass, we may replace k by the (average) neutrino 
energy E of the wave packet. This effect leads to a potentially independent way of detecting 
space-time foam [TO] , although we find that it is not competitive with limits coming directly 
from time-of-flight measurements, as discussed below. 



2.1.3 Wave-Packet Spread induced by Stochastic Fluctuations in Neutrino Ve- 
locities 

Another phenomenon that may also be induced by quantum-gravity foam models is stochastic 
fluctuation in the velocities of different neutrinos with the same energy [20J. As an example 
how this type of effect might arise from space-time foam, we consider the possibility of 
light-cone fluctuations. 

In the string-inspired models of space-time foam that we consider here, these may be 
induced by the summation over world-sheet surfaces with higher-genus topologies. These 
result in an effective stochastic fluctuation of the light cone of order [20] 



5c ~ 8 ^v^ , (6) 
M s c z 

where M s is the string scale, g s is the string coupling and E is the average energy of the 
massless (or, in the case of the neutrino, almost massless) probe. 

Light-cone fluctuations of the form ([6]) would also lead to a spread in the Gaussian wave 
packet, which is distinct from the spread induced by the refractive index (J3]). The spread 
induced by light-cone fluctuations would be linear in the quantum gravity scale, as seen 
from ([6]). Such an effect would lead to a stochastic spread in the arrival times of photons or 
neutrinos of order 

—E , A = — ^ 

cA ' 8g* 

where L is the distance of the observer from the source El We emphasize that, in contrast to 
the variation ([1]) in the refractive index - which refers to photons of different energies - the 
fluctuation (J7|) characterizes the statistical spread in the velocities of particles of the same 
energy. We note that the stochastic effect © is suppressed compared to the linear (n = 1) 
refractive index effect <§Q) by an extra power of the string coupling g s (we recall that, in the 
string model, M^lvi ck M s /g s ). 

The light-cone fluctuation effects may be thought of as inducing a time-independent 
spread a in a neutrino wave-packet that can be parametrized as (see later) 



^* = -£, A = ^r, (7) 



V(t) 

where 



(t-tn) 2 



a2 = a o + 4 1 ^— , Z = 0orl. (9) 

vlNl 



In our analysis below we may ignore effects associated with the expansion of the Universe, as we are 
dealing with neutrinos from galactic supernovae. 



5 



The expression (|Sj) is distinct from the relativistic Gaussian wave packet used above and 
in [20]. In the latter case, the stochastic light cone fluctuations already affect the spread 
|Axo| 2 at t — 0, which then may be identified with the a 2 in ([9]). 

2.2 Two-Dimensional Simulation of a Core-Collapse Supernova 

The analysis presented in this paper is based on the results of a two-dimensional, i.e., anti- 
symmetric, core-collapse simulation for a 15 M & star computed with the high-density equa- 
tion of state of Lattimer & Swesty [22]. The properties of the neutrino signal calculated in 
this model were discussed in detail in [5] . 

Unlike one- dimensional, i.e., spherically-symmetric, models, the neutrino emission during 
the post-bounce accretion phase exhibits rapid time-variability because of anisotropic mass 
flows in the accretion layer around the newly-formed neutron star. These flows are a conse- 
quence of convective overturn as well as the standing accretion-shock instability (SASI; [23]). 
which lead to large-scale, non-radial mass motions in the layer between the proto-neutron 
star surface and the accretion shock. Locally-enhanced mass infall to the compact remnant 
and asymmetric compression create hot spots that can produce transiently neutrino radi- 
ation that is more luminous and with a harder spectrum, emitted in preferred directions. 
Temporal variations of the luminosities and mean energies are expected to persist during 
the whole accretion phase, which can last hundreds of milliseconds. For electron neutrinos 
and antineutrinos, such variations could yield fractional changes of 10% and even higher 
during the most violent phases of core activity in two-dimensional models with no or only 
slow rotation [HE]. The corresponding effects are somewhat damped for muon and tau neu- 
trinos, because a smaller fraction of these neutrinos is produced in the outer layers of the 
proto-neutron star where asymmetric accretion causes the largest perturbations. 

The fluctuating neutrino emission has been shown to trigger a clearly detectable signature 
in the response of the IceCube detector in the case of a neutrino burst from a future Galactic 
supernova. Peaks of the power spectrum of the event rate are expected at the typical 
frequencies of the SASI and convective activity (between several tens of Hz and roughly 
200 Hz) [5]. 

While a softer nuclear equation of state (allowing for a more compact proto-neutron 
star) seems to favor larger signal amplitudes [5], several other sources of uncertainty in 
the model predictions need to be mentioned. One concerns the neutrino transport descrip- 
tion that is used, which even in the most sophisticated current multi-dimensional (two- or 
three-dimensional) models cannot be handled without approximations. In Refs. [5] , neutrino 
predictions were obtained in a scheme in which the full energy dependence of the transport 
problem was accounted for (including Doppler shifting, gravitational redshifting, and neu- 
trino redistribution in energy space by scattering reactions) but the two-dimensionality of 
the transport was treated in a "ray-by-ray" approximation. This means that the spatial 
transport was described by N spherically symmetric (radial) problems with N being the 
number of lateral grid zones, and it implies that the directional averaging or smoothing of 
the neutrino emission due to radiation received by an observer from different areas of the 
emitting surface is underestimated. However, in Refs. [5] also data averaged over hemi- 
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spheres (northern, southern, equatorial) were considered, which still led to easily-detectable 
signatures. Moreover, the basic effects discussed in [5] were confirmed by true multi-angle 2D 
transport results in Refs. [4] , though with the limitation of not including energy-bin coupling 
due to the effects mentioned above. We note also that flavour oscillations between v e and 
muon or tau anti-neutrinos, whose emission variations have lower amplitudes, reduces the 
modulations of the event rate in the IceCube detector only moderately. 

The largest uncertainty in current model predictions results from the two-dimensional 
nature of the most detailed simulations. Three-dimensional modelling of stellar core collapse 
is still in its infancy, and well-resolved simulations with energy- dependent neutrino transport 
are currently not available. First steps in this direction were reported in Refs. [7], but the 
models are either not evolved for interestingly long post-bounce times or the employed nu- 
merical resolution is poor and the neutrino data are not conclusive with respect to the effects 
discussed here. In constrast, Ref. [8] employed a simpler, energy-averaged ("grey"), ray-by- 
ray neutrino transport approximation and could follow the evolution of collapsing stellar 
cores in 3D over several hundred milliseconds of post-bounce accretion, through explosion, 
into the subsequent neutrino-cooling of the nascent neutron star. The radiated neutrino 
signal as visible by a distant observer was evaluated in [8] and also revealed variations with 
time on the scale of a few milliseconds, however with an amplitude of several percent only, 
instead of the 10% or more found in two-dimensional simulations. 

An analysis of the detectable consequences is in progress. Whilst in a two-dimensional 
simulation the existence of a symmetry axis directs the SASI sloshing motions of the shock 
and of the accretion flows, these motions are similar in all directions in three dimensions and 
thus appear to develop smaller amplitude in any particular direction, leading to a reduced 
fractional fluctuation of the observable neutrino emission. It should be noted, however, that 
the existing three-dimensional models still contain severe approximations and do not explore 
the range of interesting possibilities. In particular, they do not include the effects of rotation 
in the stellar core, which even for slow rates could significantly influence the growth of SASI 
spiral modes [23]. 

We base our analysis here on the more mature two-dimensional simulations, though 
adopting a somewhat optimistic point of view, in that we determine the maximal effects 
that can be expected within the detailed two-dimensional models currently available. We 
therefore consider (north-)polar emission (i.e., no averaging of luminosities and spectra over 
a wider range of latitudes) of electron antineutrinos, as predicted by the 15 M Q simulation 
with the relatively soft equation of state of Lattimer and Swesty in Refs. [5]. Possible flavor 
conversions between electron antineutrinos and heavy-lepton antineutrinos are ignored. 

2.3 Detector response to the time-varying SN neutrino signal 

As was discussed in detail in [5], IceCube or a future megaton-class water Cerenkov detector 
would be very promising for detecting the time- varying neutrino signal from a future galactic 
SN. Such detectors are designed to detect a large number of Cerenkov photons produced by 
neutrino events, and a single photon produced by a given neutrino can tag its arrival time. 
Hence the term 'event' can be used interchangeably to refer to photon or neutrino detection. 
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In the case of a SN at the fiducial distance of 10 kpc assumed here, the photon detection rate 
can be as high as ~ 10 3 /ms. This is similar to the intrinsic background rate estimated for 
IceCube. A megaton-scale water Cerenkov detector would achieve neutrino detection rates 
similar to IceCube and, in addition, would provide event-by-event information. Therefore, 
IceCube may serve as a benchmark detector for estimating a typical detection rate achievable 
for a time- varying SN neutrino signal 0. 

A schematic model of the IceCube detector response to a SN neutrino signal was used 
in [5] to estimate the detection rate, including efficiencies, for Cerenkov photons originating 
from the dominant inverse-beta reaction v e + p — > n + e + : 

n in -l L »e /10kpc\ 2 / E rms \ 2 
R - = 114 mS 10" erg s- (— ) (iFlfeV ) • (10) 
where Lp e and D are the SN luminosity and distance, respectively. The definition j5] 

-^rms = E (11) 

is used, where the average is to be taken over the neutrino distribution function. This 
estimate for the photon count rate uses an approximate inverse beta cross section of a = 
9.52 ■ 10" 44 cm 2 (£ Pe /MeV 2 ). 

We assume in our analysis that the neutrino data collected from a supernova explosion 
will consist of a list of individual neutrino events with measured energies E% and arrival 
times ti, as motivated by the fact that a low-energy water Cerenkov or scintillator detector 
is able to register the time ti of every event with high precision. The results of the simulation 
performed along the line of subsection 12.21 and described extensively in [5] are presented as a 
set of primary energy fluxes within time periods of durations ~ 3 — 5 ms, and each individual 
flux can be represented by a black-body spectrum with a given value of the mean energy. 
The fluxes are mapped into the photon counting rates using the benchmark detector response 
rate fllOl) . Knowing the mean and total energy of neutrinos leading to the photon counting 
rate in each time period, we assign statistically to each event a specific time of emission and 
energy. The distribution of one implementation of such neutrino time-energy assignments 
folded with the benchmark detector response ( flOl) is presented in Fig. [2j In order to obtain 
robust estimates of the sensitivities to novel effects in neutrino propagation as discussed 
below, we make a number of different implementations of the neutrino emission, all with 
independent statistical realizations of the thermal spectra. We then apply a technique based 
on wavelet transforms to these implementations, and analyze statistically the sensitivities to 
new physics that they may have. 



2.4 Wavelet Analysis Technique 

We use a wavelet transform technique (see [22] for a review) to analyze the neutrino time 
series generated by the simulated supernova explosion, as it is well adapted to capturing 
possible signatures of non-stationary power at many different frequencies. 

4 For the Super-Kamiokande detector with fiducial volume 22.5 kt, the corresponding neutrino detection 
rate is approximately 2 orders of magnitude smaller, but essentially background free. 
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Consider a time series, x n , n = 0, ... ,N — 1, in bins of equal width St. The wavelet 
transform is based on a wavelet function, ipo{v)i that depends on a dimensionless "time" 
parameter 77. To be admissible as a wavelet transform, this function must have zero mean 
and be localized in both time and frequency space. In choosing the wavelet function, there 
are several factors which should be considered: 

Non-orthogonality: The term "wavelet function" may be applied generically to either 
orthogonal or non-orthogonal wavelets. In orthogonal wavelet analysis, the number of con- 
volutions at each scale is proportional to the width of the wavelet basis at that scale. This 
produces a wavelet spectrum that contains discrete blocks of wavelet power, and is useful for 
signal processing as it gives the most compact representation of the signal. Unfortunately 
for time series analysis, an aperiodic shift in the time series produces a different wavelet 
spectrum. Conversely, a non-orthogonal analysis is highly redundant at large scales, where 
the wavelet spectrum at adjacent times is highly correlated. The term "wavelet basis" refers 
only to an orthogonal set of functions, and an orthogonal basis implies the use of the discrete 
wavelet transform, whereas non-orthogonal wavelet functions can be used with either discrete 
or continuous wavelet transforms. A non-orthogonal transform is useful for the analysis of 
time series where smooth, continuous variations in wavelet amplitude are expected, and is 
used in this study. 

Complexity: A real wavelet function provides only a single component, and can be used 
to isolate peaks or discontinuities. On the other hand, a complex wavelet function pro- 
vides information about both the amplitude and phase, and is better adapted for capturing 
oscillatory behavior. This is the choice made in this paper. 

Width: For concreteness, the width of a wavelet function is defined as the e-folding 
"time" of the wavelet amplitude. The resolution of a wavelet function is determined by 
the balance between the width in real space and the width in Fourier space. A narrow (in 
"time" ) function will have good time resolution but poor frequency resolution, while a broad 
function will have poor time resolution, but good frequency resolution. 

Shape: The wavelet function should reflect the type of feature present in the time series. 
For time series with sharp jumps or steps, one would choose a boxcar-like function, whereas 
for a smoothly-varying time series one would choose a smooth function such as a damped 
cosine. If one is primarily interested in wavelet power spectra, then the choice of wavelet 
function is not critical, and one function will give qualitatively similar results to another. 

Among common non-orthogonal wavelet functions, the Morlet wavelet is complex and 
contains a number of oscillations sufficient to detect narrow features of the power spectrum, 
and is the choice made here. It consists of a plane wave modulated by a Gaussian function 
in a variable r\: 

^(n) = n' l ^e^e^ 2/ \ (12) 

where uo is a dimensionless frequency. 

The continuous wavelet transform of a discrete sequence x n is defined as the convolution 
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of x n with a scaled and translated version of ipo(v) : 



(13) 



By varying the wavelet scale s and translating along the localized time index n, one can 
construct a picture showing both the amplitude of any features versus the scale and how this 
amplitude varies with time. Although it is possible to calculate the wavelet transform using 
( JT31) . it is convenient and faster to perform the calculations in Fourier space. 

To approximate the continuous wavelet transform, the convolution ( {TBI should be per- 
formed N times for each scale, where N is the number of points in the time series. By 
choosing N points, the convolution theorem allows us to perform all N convolutions simul- 
taneously in Fourier space using a discrete Fourier transform (DFT): 

xh = is&W-*^, (14) 

where k — 0, . . . , N — 1 is the frequency index. In the continuous limit, the Fourier transform 
of a function ip(t/s) is given by ip(su). According to the convolution theorem, the wavelet 
transform is the Fourier transform of the product: 

W n (s) = E™x k P{scj k )e ta "- nSt , (15) 

where u k = +fr|r and — f?j| for k < — and k > y, respectively. 

As already mentioned, following the criteria for selecting wavelets for a particular task 
described above, in this paper we process the neutrino signal using Morlet wavelets f|T2|) of 
frequency ujq, which takes the form 

Msoo) = ir-^H^e-^-^ 2 / 2 (16) 

after the Fourier transform, where H(u) is the Heaviside function: H(u) — 1 if u > 0, 
and zero otherwise. The width of this wavelet, defined as the e-folding time of the wavelet 
amplitude, is t s = v2s. The function %jj is normalized to unity: 

|&(a/)|W = l (17) 

and, in order to ensure that the wavelet transforms (1151) at all scales s are directly comparable 
to each other and to the transforms of other time series, the wavelet functions ipo at other 
scales are normalized to have unit energy: 



ip(su k ) = J ^^^(soJk). (18) 
This implies that 

VS£\4>(8U k )\ = N. (19) 



D The subscript on ip has been dropped, in order to indicate that ip has also been normalized (see later). 
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Using the convolution formula (fTBl) . the normalization of the function if> is 



{n' — n)St 


= [E^ 


(n' — n) St 


s 


V s 


s 



(20) 



and the wavelet power spectrum is defined by iH^s)) 2 . It is desirable to find a common 
normalization for the wavelet spectrum. Using the normalization in ffTB"]) . and refering to 
f fT5]) . the expectation value for | ^^(s) | 2 is equal to iV times the expectation value for \xk\ 2 . 
For a white- noise time series, this expectation value is <J 2 /N, where a 2 is the variance. Thus, 
for a white- noise process, the expectation value for the wavelet transform is |W n (s)| 2 = a 2 
for all n and s. 

Once the wavelet function is chosen, it is necessary to choose a set of scales s to use in 
(TT5]) . For our purposes, it is convenient to choose discrete scales related by powers of two: 



s j 



2-? & j 



s 0; 



J = 0,1, 



J. 



J 



Sj 



log 2 



NSt 
so 



(21) 



where so is the smallest resolvable scale and J determines the largest scale. The choice of 
a sufficiently small Sj depends on the width in spectral-space of the wavelet function. In 
the case of the Morlet wavelet, Sj ~ 0.5 is the largest value that still gives an adequate 
sampling scale. In the middle panel of Fig. [T]we use: N = 1024, St = 1.785 • 10~ 4 s, so = 2St, 
Sj = 0.125 and J = 48. We display in this panel the "cone of influence", which is indicated 
by the concave solid lines at the edges of the support of the signal: this is the region of the 
wavelet spectrum where edge effects become important, defined as the e-folding time for the 
wavelet autocorrelation power at each scale, and also gives a measure of the decorrelation 
time for a single spike in the time series. 

Since the wavelet transform is a bandpass filter with a known response function (the 
wavelet function), it is possible to reconstruct the original time series using either deconvo- 
lution or the inverse filter. In the case considered here, the reconstructed time series can be 
represented as the sum of the real parts of the wavelet transforms over all scales [26] : 



Sj 



CVo(0) 



St j RejW, 



(22) 



The factor ipo{fy removes the energy scaling while the ^J~s~j converts the wavelet transform 
to an energy density. The factor C comes from the reconstruction of a S function from its 
wavelet transform using the function ipo(i]) f[T2T) and is a constant for each type of wavelet 
function. The total energy is conserved by the wavelet transform, and the equivalent of 
Parseval's theorem for wavelet analyses is 



2 _ SjSt JV _ 1 j \W n (sj) 

° -CN^n=0^=0 — 



(23) 



To determine significance levels for either Fourier or wavelet spectra, one first needs to 
choose an appropriate background spectrum. It is then assumed that different realizations 
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of the neutrino emission process will be randomly distributed about this mean or expected 
background, and the actual spectrum is compared with this random distribution. For our 
phenomena, an appropriate background spectrum could be either white noise (with a flat 
Fourier spectrum) or red noise (increasing power with decreasing frequency). Here, for 
simplicity we choose a Gaussian white-noise background spectrum. 

We define as follows the null hypothesis for the wavelet power spectrum: we assume that 
the time series has a mean power spectrum, given simply by = 1 in case of the white 
noise. If a peak in the wavelet power spectrum appears significantly above this background 
spectrum, then it is considered to be a true feature with a certain percentage confidence. If 
x n is a normally-distributed random variable, then both the real and imaginary part of x k 
are normally distributed. Since the square of a normally distributed variable is chi-squared 
distributed with one degree of freedom (DOF), then the \xk\ 2 variable has a chi-squared 
distribution with two DOFs, \2 [2Z]- m the case that the original Fourier components are 
normally distributed, the wavelet coefficients should also be normally distributed, while the 
wavelet power spectrum | ^^.(s) | 2 should have a x| distribution. Thus, if the background 
were truly white-noise, the distribution shown in the middle panel of Fig. [1] would have a xl 
distribution for each point of (t, s). In summary, assuming a mean background spectrum of 
white (red) noise form, the distribution of the Fourier power spectrum reads: 



N\x k \ 
a 2 



PkXi (24) 



at each Fourier frequency index k, with P k being the mean spectrum value corresponding 
to the wavelet scale s at this index. (Here the sign =>- means "is distributed as".) The 
corresponding distribution for the local wavelet power spectrum is: 

\W n (s)\ 2 1 2 

— ~ 2 2 Pfe *2 (25) 

at each time n and scale s. Disregarding the relation between k and s, the relation (125]) 
is independent of the wavelet function. After finding an appropriate background spectrum, 
simply the white nose in our case, and choosing a particular confidence level for \ 2 such as 
95%, one can then use ( 1251) at each scale and build 95% contour lines, as seen in the middle 
panel of Fig. [TJ 

In order to examine fluctuations in power over a range of scales (s±, S2) (a band), one can 
define the scale-averaged wavelet power as the weighted sum of the wavelet power spectrum 
over the scales Si to s 2 : 

W 2 = SJ5t j2 \W n ( 8j )\* 
C 3=31 Sj 

The scale-averaged wavelet power can be used to examine modulation of one time series by 
another, or modulation of one frequency by another within the same time series. 

It is convenient to normalize the wavelet power by the expectation value for a white-noise 
time series. From f l26l) . this expectation value is (5j 5ta 2 )/(C S avg ), where S avg is defined as 

S avg =fe,-1 \ (27) 
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Using the normalization factor for white noise, the distribution can be modelled in analogy 
with ( 125]) . namely 

^4r\ W n ^ p 4> (28) 

where the scale-averaged theoretical spectrum is now given by 



P = S^f =n ^. (29) 



and X& is the chi-squared distribution with the number of DOFs Note that for white noise 
the normalization is such that this spectrum is still unity The number of DOFs in ( 1281) is 
modelled as 

where S m ^ = s 2°' 5 ^ 1+:i2 ^ and n a = j 2 — + The factor S avg / S mi( i corrects for the loss 
of DOFs that arises from dividing the wavelet power spectrum by the scale in (I26|) . 



3 Results of the Analysis 

3.1 Wavelet transforms of the neutrino time series 

The neutrino time series found in [51, summing over all the produced neutrino energies, 
is shown in the top panel of Fig. [1] We see that it exhibits structures on time scales 
below a hundredth of a second that appear, prima facie, to be far beyond the magnitude of 
fluctuations that could be expected from a 'featureless' white-noise spectrum. As discussed 
in the previous Section, the wavelet technique is very suitable for extracting such structures, 
and has been applied in analogous analyses of time structures in photon emissions from 
gamma-ray bursters [12] . The middle panel of Fig. [T] shows the normalized wavelet power 
spectrum, |W n (s)| 2 /<7 2 , for the time series of the neutrino emission shown in the top panel. 
The normalization by 1/cr 2 gives a measure of the power relative to white noise, and the 
colours represent the significance of the feature compared to a white-noise spectrum. We 
see that the wavelet transform picks out structures in the time series on time scales down to 
~ 2 x 10" 3 s. Several of these structures in the time series have high significance, well above 
the 95% CL for a white-noise spectrum (indicated by red contours). These can be seen in the 
bottom panel of Fig. (TJ where stuctures with time scales between 2 ms and 3 ms are selected. 
We note, in particular, the series of structures appearing at times between 0.22 and 0.34 s 
after the start. There are also structures in the band between 7 ms and 15 ms, and also 
appreciable power at longer periods. Since we are interested in obtaining the best resolution 

6 The binning is choosen at the level of a good fraction of ms, which seems to be reasonable for water 
Cerenkov and scintillator detectors. 
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possible, in the following we focus on the band corresponding to the smallest range of the 
scales where significant power is seen, namely that between 2 ms and 3 ms 0. 

Here, we investigate how these structures would be smeared out by the energy-dependent 
refractive index or by a stochastic spread in the velocities of different neutrinos with the 
same energy, the two possibilities described in the previous section. Specifically, we study 
the following possible energy dependences of the neutrino group velocity v gu : 

/ _ \ 1/2 

^ = 1±( — ] , (31) 
— = 1 ± 77^, ( 32 ) 



M, 



uLVl 



L vYM2 



2 



We also investigate a possible stochastic effect which may change the arrival times, t, of 
individual neutrinos, assuming a Gaussian probability distribution function: 



7>(t*och) = _^ exp 



1 j (34) 



2<x 2 

as discussed above. Here a = jiE 1 with constants 7; and / = 0,1, and 2. 



3.2 Sensitivity to an energy-dependent refractive index for neu- 
trinos 



Lower limits on M„ LV i, M^lvi; ^v\N2 and 7; may be calculated by requiring that the fine- 
scale time structures in the wavelet power spectrum do not disappear below the 95% CL of 
significance for a signal above the white-noise power spectrum. Specifically, for the models 
(}3~T]) . (|32|) and (1331) we apply to every neutrino event an energy-dependent time shift 

At = r t E l , (35) 

where 

Tl = d£-> (36) 

and I = 5, 1, 2. We then vary ti (M^lv?) and follow the evolution of the signal in the neutrino 
time series. If there is a non-trivial dispersive effect during propagation from the source, it 



7 However, we do note that these would be the most endangered by directional averaging of neutrinos 
streaming in all directions, as may occur in more complete 2D or 3D treatments of the neutrino transport. 
In contrast, the time structures connected to the typical SASI and convective timescales (tens of Hertz up to 
about 200 Hz) are likely to survive even in 3D (though possibly with a somewhat reduced amplitude) [71IH] ■ 
If the minimum period of neutrino variability were to increase, the sensitivity of our analysis would decrease 
correspondingly. 
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can be compensated by choosing the "correct" value of the time shift ti, in which case the 
original time structure at the source is recovered. On the other hand, dispersion at the 
source itself could not, in general, be compensated by any choice of t\. Quantitatively, the 
time structure of the supernova signal is recovered by maximizing the fraction of the scale- 
averaged power spectrum above the 95% CL line. In order to calculate a lower limit on t\ in 
any specific model, we examine the fine-scale time structures that appear above the 95% CL 
in the bottom panel of Fig. [1] and find the value of the time-shift parameter (136]) at which 
the signal above the 95% CL disappears. 

We first study a linearly energy- dependent neutrino refractive index of the form 1 + 
(E/Mvjxi). Fig. [3] displays the result of one simulation of the effect of such an energy- 
dependent refractive index, sampled in 21 bins corresponding to different time shifts T\. 
The vertical axis shows the strengths of the emissions in the structures with time scales 
between 2 x 10~ 3 s and 3 x 10 -3 s, applying a linear energy- dependent time shift T\ = 
4.2 x 10~ 5 (s/MeV). Looking at the structures that occur between 0.22 and 0.34 s after the 
start, we find that the significant portions of these small-scale structures (those that rise 
above the 95% fluctuation level for white-noise background neutrino emission) disappear 
for time delays r > 4.20 x 10~ 5 (s/MeV), corresponding to M vLV i > 2.45 x 10 13 GeV if a 
supernova distance L of 10 kpc is assumed. This sensitivity is two orders of magnitude more 
sensitive than that found in [15], namely M^lvi > 2.2 x 10 11 GeV, based on a one- dimensional 
simulation of a core-collapse supernova that did not exhibit the small time-scale structures 
seen in Fig. [U 

We have repeated this exercise with 25 different statistical realizations of the neutrino 
emission, calculating in each case the amount S of the total signal above 95% CL for different 
values of T\ sampled in 21 bins. The results of these 25 realizations can be fit quite well by a 
Gaussian distribution, as seen in Fig. [H (In this and subsequent figures, we concentrate on 
the structures with time scales between 2 x 10~ 3 s and 3 x 10~ 3 s that occur between 0.22 
and 0.34 s after the start.) One can see that the position of the maximum, which defines 
the value of r that maximizes the time structures in the signal and is expected to be zero, is 
indeed consistent with zero to within a precision of 10~ 6 (s/MeV), while the structures are 
washed out to below the lex level at 



where the number in square brackets [...] is obtained from a similar analysis of the super- 
luminal case. On the basis of this analysis, if significant time structures of the type found 
in the two-dimensional simulation [S] were to be seen in IceCube in neutrino data from a 
core-collapse supernova at a distance of 10 kpc, one could conclude that 



for a neutrino refractive index of the form 1 ± (E/M^yi). On the other hand, if no such 
structures were seen, inferring an upper limit on M^lvi would require strong independent 
confirmation of the structures found in [5], in particular by full three-dimensional simulations. 

We have repeated this analysis for the case of a quadratic energy dependence in the 
refractive index of the form 1 ± (E / M u uv2) 2 , finding the results shown in Fig. [51 In this case, 



n = 3.85 [3.95] x 10~ 5 s/MeV. 



(37) 



M^lvi > 2.68 [2.61] x 10 13 GeV. 



(38) 
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we find that the structures are washed out to below the la level when 

r 2 = 1.10 [1.11] x 1(T 6 s/MeV 2 , (39) 

where the number in square brackets is again obtained from a similar analysis of the super- 
luminal case. Hence, observation of significant time structures [5] in IceCube would imply 
that 

M„ LV2 > 0.97 [0.96] x 10 6 GeV, (40) 

if such structures were to be observed in a supernova explosion at 10 kpc, again with the 
proviso that inferring an upper limit on M^lv2 would require strong confirmation of the 
structures found in [5], specifically by full three-dimensional simulations. 

Finally, repeating this analysis for the case 1 ± E /M uLV i , we find the results shown in 
Fig. El 

n/2 = 3.10 [3.15] x 10" 4 s/VMeV, (41) 
(again, square brackets denote the superluminal case) corresponding to a sensitivity to 

M uLV i > 1.11 [1.07] x 10 22 GeV (42) 

if such structures were to be observed in a supernova explosion at 10 kpc. 



3.3 Sensitivity to a stochastic spread in neutrino velocities 

The possibility of a stochastic spread in the velocities of individual neutrinos with the same 
energy can be investigated in a similar way. As discussed in Section 2, the amount of velocity 
spread a might be energy- independent, and we consider this possibility as well as possible 
linear and quadratic energy dependences of a. As in the case of an energy- dependent time 
delay, any of these possibilities would tend to spread out and reduce the significances of the 
peaks found in the wavelet analysis shown in Fig. [TJ 

This effect is seen for a oc E in Fig. and for a oc E 2 in Fig. [9j In each case, we plot 
results obtained from 35 independent statistical simulations of the neutrino emission signal. 
We see that the wavelet peaks are reduced below the 95% CL white-noise level for 

r stochl = 2.16 x 10- 5 s/MeV (43) 

in the linear case, 

r s toch 2 = 9.56 x KT 7 s/MeV 2 (44) 

in the quadratic case, and 

r stoch0 = 3.59 x 10- 4 s (45) 

in the energy- independent caseB These sensitivities correspond to M u ^ 1 > 4.78 x 10 13 GeV 
and M v ^f 2 > 1.04 x 10 6 GeV in the energy- dependent cases, for a supernova explosion at 
10 kpc. 

8 For comparison, we note that the OPERA result [17] also provides an upper limit on the stochastic 
spread of about 10 ns, though after propagation over a much shorter distance ~ 730 km. However, using 
equations ([2} to ([5|) we find that neither OPERA nor SN1987a gives competitive constraints on the scales 
M vVn . 
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4 Conclusions and Prospects 



We have shown that the existence of structures with short time scales in the neutrino emis- 
sion from a core-collapse supernova, as suggested by two-dimensional simulations [5], would 
open up new possibilities for probing aspects of the propagation of neutrinos that lie far 
beyond the reach of terrestrial experiments, and up to two orders of magnitude beyond the 
sensitivity provided by previous analyses based on one- dimensional supernova simulations. 
This increased sensitivity holds for possible square-root, linear and quadratic dependences of 
the neutrino refractive index, and for both energy-independent and linear or quadrat ically- 
dependent stochastic spreads in the velocities of different neutrinos with the same energy. 

Specifically, if such short time structures are seen in a supernova explosion at a distance 
of 10 kpc, one could infer that 



M vLv i > 1.11 [1.07] x 10 22 GeV, (46) 

ilWi > 2.68 [2.61] x 10 13 GeV, (47) 

M uLV2 > 0.97 [0.96] x 10 6 GeV, (48) 

M vWl > 4.78 x 10 13 GeV, (49) 

M vm > 1.04 x 10 6 GeV, (50) 



where the numbers in square brackets correspond to the superluminal case, and the last 
two limits correspond to the possible effects of stochastic fluctuations. In the case of an 
energy-independent stochastic spread, one could infer that r stoc ho < 3.59 x 10™ 4 s. 

If such short time structures are not seen, many checks would be necessary before one 
could conceivably claim observation of any unconventional effect in neutrino propagation. 
In particular, it would be necessary to validate the predictions of the two-dimensional core- 
collapse supernova simulation on which this analysis is based, specifically by confirming that 
short time structures are also found in full three-dimensional simulations [7J[S]. We hope that 
the interesting sensitivity to new neutrino physics discussed in this paper - not to mention 
OPERA [17] - will add to the motivation to develop further such simulations and derive 
robust predictions for neutrino emissions from core-collapse supernovae. 
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Figure 1: Top panel: The time series of the neutrino emission from the two-dimensional 
simulation of a core-collapse supernova found in |5]/. The time profile is sampled in 1024 
(2 10 ) bins. Middle panel: The local wavelet power spectrum of the neutrino emission time 
series, obtained using the Morlet wavelet function normalized by 1/a 2 . The vertical 
axis is the Fourier period (in seconds), and the horizontal axis is the time of the neutrino 
emission. The red contours enclose regions that differ from white noise at greater than the 
95% confidence level. The cone of influence, where edge effects become important, is indicated 
by the concave solid lines at the edges of the support of the signal. Comparing the width of a 
peak in the wavelet power spectrum with this decorrelation time, one can distinguish between a 
spike in the data (possibly due to random noise) and a harmonic component at the equivalent 
Fourier frequency. Bottom panel: The average power in the 0.002 - 0.003 s band. The dashed 
line is the 95% confidence level obtained from K2~tyl . 
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Figure 2: The distribution of times and energies assigned to individual neutrinos in one 
statistical realization of the thermal spectra found in the simulation J^j. 
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Figure 3: The strength of the time-scale structure of the power spectrum averaged between 
0.002 and 0.003 s disappears below the 95% CL of significance after applying a linear energy- 
dependent time shift T\ = 4.2 x 1CT 5 (s/MeV). 
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Figure 4: A Gaussian fit to the amount E of the short time-scale signal above the 95% CL, 
calculated for 21 values of the shift parameters T\. Each point is obtained as the average over 
25 realizations of the time-energy assignments of individual neutrinos. 
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Figure 5: A Gaussian fit to the amount E of the short time-scale signal above the 95% CL 
calculated for 21 values of the shift parameters T2- Each point is obtained as the average over 
25 realizations of the time-energy assignments of individual neutrinos. 
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Figure 6: A Gaussian fit to the amount E of the short time-scale signal above the 95% CL 
calculated for 26 values of the shift parameters ri . Each point is obtained as the average 
over 25 realizations of the time-energy assignments of individual neutrinos. 
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Figure 7: A Gaussian fit to the amount E of the short time-scale signal above the 95% CL 
calculated for 11 values of the shift parameters rf toch . Each point is obtained as the average 
over 35 realizations of the time-energy assignments of individual neutrinos. 
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Figure 8: A Gaussian fit to the amount E of the short time-scale signal above the 95% CL 
calculated for 17 values of the shift parameters T"o toch . Each point is obtained as the average 
over 35 realizations of the time-energy assignments of individual neutrinos. 
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Figure 9: A Gaussian fit to the amount E of the short time-scale signal above the 95% CL 
calculated for 13 values of the shift parameters r| toch . Each point is obtained as the average 
over 35 realizations of the time-energy assignments of individual neutrinos. 
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